packages = c('rgdal', 'maptools', 'raster','spatstat', 'tmap', 'sf', 'tidyverse', 'SpatialAcc', 'tbart', 'cclust', 'rgeos')
for (p in packages){
if(!require(p, character.only = T)){
install.packages(p) }
library(p,character.only = T) }
Singapore Planning Subzone (MP14_SUBZONE_WEB_PL)
mpsz = st_read(dsn = "data", layer = "MP14_SUBZONE_WEB_PL")
Reading layer `MP14_SUBZONE_WEB_PL' from data source `C:\IS415\IS415-Neighbourhood-WatchDocs\data' using driver `ESRI Shapefile'
Simple feature collection with 323 features and 15 fields
geometry type: MULTIPOLYGON
dimension: XY
bbox: xmin: 2667.538 ymin: 15748.72 xmax: 56396.44 ymax: 50256.33
epsg (SRID): NA
proj4string: +proj=tmerc +lat_0=1.366666666666667 +lon_0=103.8333333333333 +k=1 +x_0=28001.642 +y_0=38744.572 +datum=WGS84 +units=m +no_defs
mpsz <- mpsz[order(mpsz$SUBZONE_N),]
mpsz1 <- readOGR(dsn = "data", layer = "MP14_SUBZONE_WEB_PL")
OGR data source with driver: ESRI Shapefile
Source: "C:\IS415\IS415-Neighbourhood-WatchDocs\data", layer: "MP14_SUBZONE_WEB_PL"
with 323 features
It has 15 fields
mpsz1
class : SpatialPolygonsDataFrame
features : 323
extent : 2667.538, 56396.44, 15748.72, 50256.33 (xmin, xmax, ymin, ymax)
coord. ref. : +proj=tmerc +lat_0=1.366666666666667 +lon_0=103.8333333333333 +k=1 +x_0=28001.642 +y_0=38744.572 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0
variables : 15
names : OBJECTID, SUBZONE_NO, SUBZONE_N, SUBZONE_C, CA_IND, PLN_AREA_N, PLN_AREA_C, REGION_N, REGION_C, INC_CRC, FMEL_UPD_D, X_ADDR, Y_ADDR, SHAPE_Leng, SHAPE_Area
min values : 1, 1, ADMIRALTY, AMSZ01, N, ANG MO KIO, AM, CENTRAL REGION, CR, 00F5E30B5C9B7AD8, 2014/12/05, 5092.895, 19579.07, 871.5549, 39437.94
max values : 323, 17, YUNNAN, YSSZ09, Y, YISHUN, YS, WEST REGION, WR, FFCCF172717C2EAF, 2014/12/05, 50424.792, 49552.79, 68083.9365, 69748298.79
subzone <- mpsz[mpsz$SUBZONE_N == "ADMIRALTY",]
p <- as(subzone, 'Spatial')
sp_cent <- gCentroid(p, byid = TRUE)
#crs(sp_cent)
# sp_cent <- st_transform(subzone, crs = st_crs(sp_cent))
#
sp_cent
class : SpatialPoints
features : 1
extent : 27091.28, 27091.28, 48333.8, 48333.8 (xmin, xmax, ymin, ymax)
coord. ref. : +proj=tmerc +lat_0=1.366666666666667 +lon_0=103.8333333333333 +k=1 +x_0=28001.642 +y_0=38744.572 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0
#plot(subzone)
Singapore Residents by Subzone, Age Group and Sex (Gender)
population = read_csv("data/respopagsex2000to2017.csv")
Parsed with column specification:
cols(
PA = [31mcol_character()[39m,
SZ = [31mcol_character()[39m,
AG = [31mcol_character()[39m,
Sex = [31mcol_character()[39m,
Pop = [32mcol_double()[39m,
Time = [32mcol_double()[39m
)
Singapore Residents by Subzone, Age Group and Sex, June 2017 (Gender)
popage2017pre <- population %>%
filter(Time == 2017) %>%
filter(AG == "65_to_69" | AG == "70_to_74" | AG == "75_to_79" | AG == "80_to_84" | AG == "85_and_over") %>%
group_by(SZ)
popage2017pre_added <- aggregate(Pop ~ SZ, popage2017pre, sum)
Total Number of Clinics in Singapore
clinics = read_csv("data/clinics_hcidirectory.csv")
Parsed with column specification:
cols(
clinic_name = [31mcol_character()[39m,
address = [31mcol_character()[39m,
postal_code = [31mcol_character()[39m,
X = [32mcol_double()[39m,
Y = [32mcol_double()[39m,
LAT = [32mcol_double()[39m,
LONG = [32mcol_double()[39m
)
clinics_essentials <- c("clinic_name", "address", "LAT", "LONG", "X", "Y") # remove postal code
clinics <- clinics[clinics_essentials]
Total Number of TCM Clinics in Singapore
tcm = read_csv("data/tcm_tcmboard.csv")
Parsed with column specification:
cols(
tcm_physician_name = [31mcol_character()[39m,
registration_no = [31mcol_character()[39m,
tcm_place_name = [31mcol_character()[39m,
tcm_address = [31mcol_character()[39m,
postal_code = [31mcol_character()[39m,
X = [32mcol_double()[39m,
Y = [32mcol_double()[39m,
LAT = [32mcol_double()[39m,
LONG = [32mcol_double()[39m
)
tcm <- tcm[!duplicated(tcm$tcm_place_name), ]
tcm_essentials <- c("tcm_place_name", "tcm_address", "LAT", "LONG", "X", "Y")
tcm <- tcm[tcm_essentials]
names(tcm)[names(tcm) == "tcm_place_name"] <- "clinic_name"
names(tcm)[names(tcm) == "tcm_address"] <- "address"
Combine TCM and clinics
clinics_combined <- rbind(clinics, tcm)
clinics_combined <- st_join(st_as_sf(clinics_combined,
coords = c("X", "Y"), crs = st_crs(mpsz)), mpsz)
Number of HDB blocks per planning subzone
HDB = read_csv("data/hdb_property_information.csv")
Parsed with column specification:
cols(
blk_no_street = [31mcol_character()[39m,
total_dwelling_units = [32mcol_double()[39m,
`1room_sold` = [32mcol_double()[39m,
`2room_sold` = [32mcol_double()[39m,
`3room_sold` = [32mcol_double()[39m,
`4room_sold` = [32mcol_double()[39m,
`5room_sold` = [32mcol_double()[39m,
exec_sold = [32mcol_double()[39m,
multigen_sold = [32mcol_double()[39m,
studio_apartment_sold = [32mcol_double()[39m,
`1room_rental` = [32mcol_double()[39m,
`2room_rental` = [32mcol_double()[39m,
`3room_rental` = [32mcol_double()[39m,
other_room_rental = [32mcol_double()[39m,
postal_code = [31mcol_character()[39m,
X = [32mcol_double()[39m,
Y = [32mcol_double()[39m,
LAT = [32mcol_double()[39m,
LONG = [32mcol_double()[39m
)
HDB_sf <- st_as_sf(HDB, coords = c("X", "Y"), crs = st_crs(mpsz))
Residents by Age Group & Type of Dwelling, Annual
popByDwelling = read_csv("data/residents-by-age-group-and-type-of-dwelling-detailed-categories-annual.csv")
Parsed with column specification:
cols(
year = [32mcol_double()[39m,
level_1 = [31mcol_character()[39m,
level_2 = [31mcol_character()[39m,
level_3 = [31mcol_character()[39m,
value = [32mcol_double()[39m
)
Residents by Age Group & Type of Dwelling, Annual 2017
popByDwelling2017 <- popByDwelling %>%
filter(year == 2017) %>%
filter(level_1 == "65-69 Years"
| level_1 == "85 Years & Over"
| level_1 == "70-74 Years"
| level_1 == "75-79 Years"
| level_1 == "80-84 Years") %>%
group_by(level_3)
popByDwelling2017
popByDwelling2017_added <- aggregate(value ~ level_3, popByDwelling2017, sum)
No. of blocks per subzone
mpsz_HDB <- st_join(HDB_sf,mpsz)
mpsz_HDB$`1_2room_total` <- mpsz_HDB$`1room_sold` + mpsz_HDB$`2room_sold` + mpsz_HDB$`1room_rental` + mpsz_HDB$`2room_rental`
mpsz_HDB_1_2_room = aggregate(mpsz_HDB$`1_2room_total`, by=list(SUBZONE=mpsz_HDB$SUBZONE_N), FUN=sum) %>%
rename('No_of_1_2_room' = 'x')
mpsz_HDB$`3room_total` <- mpsz_HDB$`3room_sold` + mpsz_HDB$`3room_rental` + mpsz_HDB$`other_room_rental`
mpsz_HDB_3_room_added = aggregate(mpsz_HDB$`3room_total`, by=list(SUBZONE=mpsz_HDB$SUBZONE_N), FUN=sum) %>%
rename('No_of_3_room' = 'x')
mpsz_HDB$`4room_total` <- mpsz_HDB$`4room_sold`
mpsz_HDB_4_room_added = aggregate(mpsz_HDB$`4room_total`, by=list(SUBZONE=mpsz_HDB$SUBZONE_N), FUN=sum) %>%
rename('No_of_4_room' = 'x')
mpsz_HDB$`5room_exec_total` <- mpsz_HDB$`5room_sold` + mpsz_HDB$`exec_sold`
mpsz_HDB_5_room_exec_added = aggregate(mpsz_HDB$`5room_exec_total`, by=list(SUBZONE=mpsz_HDB$SUBZONE_N), FUN=sum) %>%
rename('No_of_5_room_exec' = 'x')
mpsz_HDB_added <- left_join(mpsz_HDB_1_2_room, mpsz_HDB_3_room_added, by ='SUBZONE') %>%
left_join(., mpsz_HDB_4_room_added, by = 'SUBZONE') %>%
left_join(., mpsz_HDB_5_room_exec_added, by = 'SUBZONE')
mpsz_HDB_added$`total_units` <- mpsz_HDB_added$`No_of_1_2_room` + mpsz_HDB_added$`No_of_3_room` + mpsz_HDB_added$`No_of_4_room` + mpsz_HDB_added$`No_of_5_room`
total_1_2_room_units <- sum(mpsz_HDB_added$No_of_1_2_room)
total_3_room_units <- sum(mpsz_HDB_added$No_of_3_room)
total_4_room_units <- sum(mpsz_HDB_added$No_of_4_room)
total_5_room_exec_units <- sum(mpsz_HDB_added$No_of_5_room_exec)
No. of elderly per block in every subzone
total_1_2_room_units <- sum(mpsz_HDB_added$No_of_1_2_room)
total_3_room_units <- sum(mpsz_HDB_added$No_of_3_room)
total_4_room_units <- sum(mpsz_HDB_added$No_of_4_room)
total_5_room_units <- sum(mpsz_HDB_added$No_of_5_room)
mpsz_HDB$No_of_Elderly_in_block_1_2 <- ifelse(mpsz_HDB$`1_2room_total` == 0, 0, mpsz_HDB$`1_2room_total`/total_1_2_room_units) * popByDwelling2017_added$value[popByDwelling2017_added$level_3=="HDB 1- And 2-Room Flats"]
mpsz_HDB$No_of_Elderly_in_block_3 <- ifelse(mpsz_HDB$`3room_total` == 0, 0, mpsz_HDB$`3room_total`/total_3_room_units) * popByDwelling2017_added$value[popByDwelling2017_added$level_3=="HDB 3-Room Flats"]
mpsz_HDB$No_of_Elderly_in_block_4 <- ifelse(mpsz_HDB$`4room_total` == 0, 0, mpsz_HDB$`4room_total`/total_4_room_units) * popByDwelling2017_added$value[popByDwelling2017_added$level_3=="HDB 4-Room Flats"]
mpsz_HDB$No_of_Elderly_in_block_5 <- ifelse(mpsz_HDB$`5room_exec_total` == 0, 0, mpsz_HDB$`5room_exec_total`/total_5_room_exec_units) * popByDwelling2017_added$value[popByDwelling2017_added$level_3=="HDB 5-Room And Executive Flats"]
mpsz_HDB$No_of_Elderly_in_block <- mpsz_HDB$No_of_Elderly_in_block_1_2 + mpsz_HDB$No_of_Elderly_in_block_3 + mpsz_HDB$No_of_Elderly_in_block_4 + mpsz_HDB$No_of_Elderly_in_block_5
mpsz_HDB$No_of_Elderly_in_block <- round(mpsz_HDB$No_of_Elderly_in_block)
#assign capacity of 2 per clinic
clinics_combined['capacity'] <- 2
mpsz_HDB_coords <- mpsz_HDB %>% st_coordinates()
clinics_combined_coords <- clinics_combined %>% st_coordinates()
dm <- distance(mpsz_HDB_coords, clinics_combined_coords)
acc_hansen <- data.frame(ac(mpsz_HDB$No_of_Elderly_in_block, clinics_combined$capacity, dm, power=0.01, family="Hansen"))
colnames(acc_hansen) <- "accHansen"
acc_hansen <- tbl_df(acc_hansen)
HDB_acc <- bind_cols(mpsz_HDB, acc_hansen)
tmap_mode("view")
tmap mode set to interactive viewing
tm_shape(mpsz_HDB)+
tm_symbols(size = 0.1)+
tm_shape(HDB_acc)+
tm_bubbles(col = "accHansen",
n=5,
style = "quantile",
size = 0.1,
border.col = "black",
border.lwd = 1)